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Abstract 

We  derive  a  semi-empirical  equation  to  describe  the  performance  curves  of  polymer  electrolyte  membrane  fuel  cells  (PEMFCs).  The 
derivation  is  based  on  the  observation  that  the  main  non-linear  contributions  to  the  cell  voltage  deterioration  of  Hy/air  feed  cells  are  deriving 
from  the  cathode  reactive  region.  To  evaluate  such  contributions  we  assumed  that  the  diffusion  region  of  the  cathode  is  made  by  a  network  of 
pores  able  to  transport  gas  and  liquid  mixtures,  while  the  reactive  region  is  made  by  a  different  network  of  pores  for  gas  transport  in  a  liquid 
permeable  matrix.  The  mathematical  model  is  largely  mechanistic,  with  most  terms  deriving  from  phenomenological  mass  transport  and 
conservation  equations.  The  only  full  empirical  term  in  the  performance  equation  is  the  Ohmic  overpotential,  which  is  assumed  to  be  linear 
with  the  cell  current  density.  The  resulting  equation  is  similar  to  other  published  performance  equations  but  with  the  advantage  of  having 
coefficients  with  a  precise  physical  origin,  and  a  precise  physical  meaning.  Our  semi-empirical  equation  is  used  to  fit  several  set  of  published 
experimental  data,  and  the  fits  showed  always  a  good  agreement  between  the  model  results  and  the  experimental  data.  The  values  of  the  fitting 
coefficients,  together  with  their  associated  physical  meaning,  allow  us  to  asses  and  quantify  the  phenomenology  which  is  set  on  in  the  cathode 
as  the  cell  current  density  is  increased.  More  precisely,  we  observe  the  development  of  the  flooding  and  of  the  local  decrease  of  the  oxygen 
concentration.  Further  developments  of  such  a  model  for  the  cathode  compartment  of  the  fuel  cell  are  discussed.  ©  2002  Elsevier  Science 
B.V.  All  rights  reserved. 
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1.  Introduction 

The  last  decade  saw  an  explosion  in  the  number  of  papers 
dealing  with  the  modelling  of  fuel  cells,  and  in  particular 
with  the  modelling  of  polymer  electrolyte  membrane  fuel 
cells  (PEMFCs).  For  what  concerns  the  spatial  dimension¬ 
ality,  published  models  include  0D  models  [  1—6],  in  which 
no  spatial  dimensions  of  the  cell  are  explicitly  taken  into 
account,  ID  models  [7-11],  in  which  the  spatial  dimension 
entering  into  the  description  is  the  one  perpendicular  to  the 
plates  of  the  cell,  and  2D  models  [12-15],  in  which  the 
considered  planes  are  those  perpendicular  to  the  cell  plates. 
To  our  knowledge,  3D  models  in  which  the  single  cell  is 
considered  in  its  full  glory  have  been  only  published  for 
solid  oxide  fuel  cells  [16],  For  what  concerns  the  temporal 
dimensionality,  both  steady-state  and  dynamical  models 
already  appeared.  While  for  what  concerns  the  included 
phenomenology,  the  larger  amount  of  published  models 
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deals  with  charge  and  mass  transport  and  conservation. 
Few  also  include  heat  transport. 

Clearly,  the  complexity  of  the  models  is  extremely  vari¬ 
able,  and  it  goes  from  the  one-equation  approach  typical  of 
the  empirical  models,  to  very  intimidating  set  of  coupled 
differential  equations  with  tens  of  spatial-dependent  mate¬ 
rial  properties.  It  is  believed  that  the  pay  back  of  complex 
models  is  a  strong  increase  in  their  predictive  power.  How¬ 
ever,  and  before  to  develop  new  complex  approaches,  it  is 
always  better  to  fully  exploit  the  potentialities  of  more 
simple  point  of  view. 

An  open  problem  of  simple  empirical  approaches  which 
try  to  mimic  the  performance  curve  of  PEMFC  is  the 
physical  origin  and  meaning  of  the  cell  current-dependent 
terms.  This  is  true  for  some  of  the  most  spread-out  perfor¬ 
mance  equations  which  are  that  of  Kim  et  al.  [1] 

kcell  =  ^O.K  —  Rcell.K/  ~  b K  In  (/)  —  mK  eXp(«K/),  (1) 

that  of  Lee  et  al.  [2] 

Vceii  =  £o,L  -  Rceii+Z  -  bL  In  (/)  -  mL  exp  (nLI) 

-aL  In  (Pj Pq2),  (2) 
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Nomenclature 

a  effective  catalyst  area  per  unit  volume 

Co,  concentration  of  oxygen 

Cd  cathode  diffusion  region 

Cr  cathode  reaction  region 

Dj-j  diffusion  coefficient  of  the  pair  i—j  in  a  gas 
binary  mixture 

Do,  diffusion  coefficient  of  oxygen 

Dp  mean  pore-pore  distance 

F  Faraday  constant 

i  ionic  current  density 

z'o  exchange  current  density 

I  cell  current  density 

K  proportionality  constant 

Kt]  Henry’s  law  constant 

L  diffusion  length  of  oxygen 

Lc  characteristic  diffusion  length 

L\  longitudinal  diffusion  length 

Lp  mean  pore  length 

Lt  transversal  diffusion  length 

M  membrane  region 

Ni  diffusion  mechanism  parameter 

No2  mass  flux  of  oxygen 

P  hydraulic  pressure 

R  gas  constant 

/?ceii  cell  resistence 

Rp  mean  pore  radius 

,S’e  number  of  electrons  exchanged  in  the  electro¬ 

chemical  reaction 
S  flooding  parameter 

T  temperature 

Vceii  cell  potential 

VcR  volume  of  Cr  region 

Xj  molar  fraction  of  species  i 

Greek  letters 

aa  anode  transfer  coefficient 

ac  cathode  transfer  coefficient 

p F  Faraday  constant  in  units  of  RT 

e„  gas  porosity 

y  kinetic  exponent  of  the  species  in  the  Butler- 

Volmer  equation 
t]  overpotential 

r/)|  electric  potential  of  the  liquid  phase 

<j)s  electric  potential  of  the  solid  phase 

p  empirical  constant 

cr  area  crossed  by  oxygen  to  make  its  electro¬ 

chemical  reaction 
crc  cross-area  of  the  cell 

E  integrated  oxygen  concentration 

Subscripts 

int  interface  value 

Ohm  Ohmic  value 


ref  at  the  reference  state 
w  water 

Superscripts 

0  value  at  the  zero-current  density 

eff  effective  in  porous  media 

1  value  at  the  limiting  current  density 


and  that  of  Squadrito  et  al.  [3] 

Feel  i  =  E0,s  -  /?CC||,S/  -  bs  In  (/)  +  asIk  In  ( 1  -pi).  (3) 

Here,  Vceii  is  the  cell  potential,  I  the  cell  current  density,  P 
the  pressure,  and  all  the  other  quantities  are  the  fitting 
coefficients. 

Since  all  these  equations  are  rich  of  fitting  coefficients, 
they  always  reproduce  experimental  data  in  an  accurate  way. 
However,  very  little  is  learned  about  the  phenomena  which 
take  place  in  the  cell  as  the  cell  current,  /,  changes.  In 
addition,  the  systematic  application  of  these  equations  does 
not  show  clear  trends  in  the  fitting  coefficients,  eliminating 
the  possibility  of  using  them  as  predictive  tools. 

To  repair  this  state  of  affair,  Amphlett  and  coworkers 
[4-6]  made  a  semi-empirical  derivation  of  a  performance 
equation,  in  which  both  mechanistic  and  empirical  features 
are  present  in  the  fitting  coefficients.  However,  the  proposed 
equation  shows  some  lacks  in  reproducing  experimental 
data,  especially  at  high  current  density. 

In  the  same  spirit,  we  also  present  a  semi-empirical 
derivation  of  a  performance  equation  with  the  goal  of  having 
the  largest  number  of  mechanistic  derived  coefficients.  The 
hope  is  that  once  the  equation  is  used  to  fit  experimental 
data,  the  values  of  the  mechanistic  derived  coefficients  will 
give  us  information  on  the  phenomena  which  are  set  on  in 
the  cell  as  the  current  is  changed. 

To  make  the  semi-empirical  derivation  of  the  performance 
equation,  we  pay  attention  to  transport  cathode  phenomena 
by  introducing  a  porous  structure  in  both  the  diffusion  and 
reactive  regions  of  the  cathode  itself.  This  choice  is  made  to 
introduce  into  the  model  flooding  and  concentration  effects, 
which  severely  affect  the  fuel  cell  behavior  at  high  cell 
current  density. 

After  an  introductory  and  qualitative  analysis  of  perfor¬ 
mance  curve  made  in  Section  2,  Section  3  of  this  paper  is 
devoted  to  the  development  of  our  model.  Section  4  shows 
comparisons  with  other  published  performance  equations, 
with  experimental  data,  and  with  some  of  the  published 
mechanistic  models,  too.  Section  5,  contains  summary  and 
conclusions. 


2.  Qualitative  analysis  of  performance  curves 

Potential  losses  inside  a  fuel  cell  originate  either  at  the 
interfaces  between  solid  and  fluid  phase,  where  all  the 
electrochemistry  takes  place,  or  internally  to  each  single 
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ad  Ar  m  cr  cd 


Fig.  1 .  Typical  cell  potential  drop  within  the  different  regions  of  a  PEMFC 
at  a  given  cell  current  density.  Ad,  Cd:  anode  and  cathode  diffusion 
regions;  Ar,  Cr:  anode  and  cathode  reaction  regions;  M:  membrane; 
4>i,(j)s:  liquid  and  solid  phase  potentials. 


phase.  The  single  phase  losses  are  essentially  of  Ohmic 
origin,  linear  with  the  cell  current  density,  7,  and  therefore, 
they  have  a  trivial  effect  on  the  performance  curves  (cell 
voltage  against  cell  current  density).  More  complex  is  the  7- 
behavior  of  the  interface  overpotentials  which,  as  illustrated 
in  Fig.  1,  originate  in  the  small  electrode  reactive  regions 
where  liquid  phase  and  solid  phase  coexist. 

Two  main  factors  govern  the  magnitude  of  the  interface 
overpotentials:  the  reaction  kinetic  and  the  local  concentra¬ 
tion  or,  in  other  words,  the  local  availability  of  reactants. 

As  appear  in  Fig.  1,  the  anodic  term  is  much  smaller 
than  the  cathodic  one.  This  is  because  of  the  extremely  favo¬ 
rable  electrochemical  kinetic  constant  associated  with  the 


hydrogen  reaction  as  compared  to  that  of  the  oxygen  reac¬ 
tion,  and  because  hydrogen  supply  do  not  present  problems. 
Indeed,  hydrogen  is  supplied  as  pure  H2  (or  with  small  con¬ 
centrations  of  contaminants),  while  oxygen  is  mostly  sup¬ 
plied  as  air.  Furthermore,  flooding  of  the  diffusive  regions  is 
more  severe  in  the  cathode  than  in  the  anode,  because  liquid 
water  is  generated  during  the  cathode  reaction. 

The  above  discussion  is  also  well  illustrated  by  Fig.  2 
which  shows  typical  results  obtained  from  simulation  meth¬ 
ods  (in  this  case,  results  are  from  the  modified  Bernardi- 
Verbrugge  model  [8,9]).  The  contributions  of  the  different 
regions  of  the  PEMFC  to  the  cell  potential,  Vceii,  are  now 
plotted  as  a  function  of  the  cell  current  density,  7.  It  is  seen 
that  the  main  non-linear  contribution  comes  from  the  cath¬ 
ode  reactive  region,  Cr,  which  accounts  for  both  the  activa¬ 
tion  and  concentration  overpotentials. 

Therefore,  we  write  the  cell  potential  as  made  up  of  a 
zero-current  term,  £’(7  =  0),  MEA  potential  drop  terms  of 
Ohmic  origin,  f/ohmCO’  and  potential  drop  terms  originated 
by  the  interface  electrochemistry  in  the  anode  and  cathode 
reactive  regions,  >/“nt(7): 

VcenW  =  E(I  =  0)  +  ?/ohm(7)  +  r,  ghm(7)  +  rfi&Jt) 

+  '/kit  CO  +  '7int(A  (4) 

and  we  restrict  our  analysis  to  the  j/gt(7)  term  by  rewriting 
Vceii(7)  as 

ycell(7)  ~  E{I  =  0)  -  /W  +  {!)■  (5) 

Extension  of  the  analysis  to  i] ^t(7)  can  be  done  by  following 
the  same  lines  with  only  minor  changes. 


Current  density  (A/cm2) 


Fig.  2.  Performance  curve  and  cell  overpotentials. 
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3.  The  semi-empirical  approach 


3.1.  Total  oxygen  concentration 


To  relate  the  interphase  cathodic  overpotential,  to 

the  cell  current  density,  7,  we  start  from  the  Butler- Volmer 
equation 

e«aft>(<Mz)-0i(z))  _  e-«<A(0s(z)-0 l(z)) 

(6) 

where  i  is  the  ionic  current  density,  a  the  effective  catalyst 
area  per  unit  of  volume,  *"o,ref  the  exchange  current  density, 
Co,  the  oxygen  concentration,  y  a  kinetic  coefficient,  a,  the 
electrode  transfer  coefficients,  and  />F  the  Faraday  constant 
in  units  of  RT. 

The  values  taken  by  (f>s  and  0,  within  the  region  Cr  (see 
Fig.  1)  indicate  that  the  first  exponential  of  Eq.  (6)  (the  anode 
semi-reaction)  is  completely  negligible,  and  that  the  poten¬ 
tial  variations  of  Ohmic  origin  are  negligible,  too,  due  to  the 
small  dimensions  of  the  region.  Within  a  good  approxima¬ 
tion,  Eq.  (6)  can  be  re-written  as 


V  •  i  =  aio 


ref 


C02(Z) 


V  •  i  «  -aio, ref 
,„c 


Co2 


=-“c  A0 


f02.rcl 

where  t]\nl  is  defined  as 


P„t  =  ^R/CD-^R/CD- 


'?Pnt 


(7) 

(8) 


The  integration  of  the  Butler- Volmer  equation  over  the 
cathode  reactive  region  gives  the  total  amount  of  current 
produced  in  the  cell,  i.e.  the  cell  current  density  I.  Then,  by 
integrating  Eq.  (7)  in  Cr  we  get 


acI=  [  V  •  ;dV  «  -  f  dy;  (9) 

-'Cr  ^02, ref  7Cr 

in  which,  crc  is  the  cross-area  of  the  cell. 

From  Eq.  (9)  and  by  throwing  all  the  constant  terms  in  a 
generic  constant  K ,  we  get  for  the  expression 


'lint  CO  =  K  —  — [  In  (/)  —  In  (£(/))],  (10) 

«cPf 

where  we  have  defined 


s(i)=  [  (C02(i)y dv.  (id 

JcR 

At  this  point,  the  problem  is  shifted  to  the  finding  of  an 
expression  for  3(7).  This  is  done  in  two  steps.  First,  3(7)  is 
expressed  as  a  function  of  the  interface  concentration 
c§fD(7).  Second,  an  expressions  for  c'q*/Cd (!)  is  found. 
The  first  step  involves  the  mechanisms  of  oxygen  diffusion 
inside  the  Cr  region,  while  the  second  step  concerns  oxygen 
diffusion  in  the  Cd  region.  Because  3(7)  in  Eq.  (10)  is  inside 
a  logarithm,  all  the  multiplying  factors  which  do  not  depend 
on  7  (including  numerical  factors,  physical  constants,  mea¬ 
sure  units,  etc.)  can  be  thrown  away  in  the  trashcan  constant 
K,  and  we  can  proceed  by  proportionality  relations. 


In  the  reactive  region,  the  material  conservation  equation 
for  oxygen  reads  as  V  •  Vo,  cx  V  •  i.  Such  an  equation  can  be 
integrated  in  Cr,  and  by  using  the  condition  of  negligible 
oxygen  flux  at  M/Cr  interface  and  the  first  equality  in 
Eq.  (9),  we  have  ct|Vo^Cd|  oc  (xc7.  By  using  this  result 
and  the  phenomenological  diffusion  equation, 
Vo,  =  —  Dq2  Vco,,  the  gradient  of  the  oxygen  concentration 
at  the  Cr/Cd  interface  becomes 

l^cO,  Icr/Cd  ^  “  •  (12) 

Here,  we  made  explicit  the  dependence  on  cr,  which  is  the 
surface  crossed  by  oxygen  to  make  its  electrochemical 
reaction,  since  it  can  change  according  to  the  physical  model 
chosen  to  describe  oxygen  diffusion  within  the  reactive 
region. 

By  assuming  a  linear  decrease  of  the  oxygen  concentra¬ 
tion  profile  within  the  Cr  region,  the  length  7,(7)  over  which 
oxygen  can  diffuse  (i.e.  the  distance  over  which  cq,(7)  0  0) 
is  defined  as 


7.(7) 


CofD(I)  c£/Cd(7) 

— — - £--7 - oc  o — - - 

|Vco,(7)|Cr/Cd  7 


(13) 


where  we  also  made  use  of  Eq.  (12). 

Now,  to  evaluate  the  volume  integral  of  the  linear  function 


«„.(/,/)  =  FS,C”(')(1-W/))  for  1<W). 

lo  for  l  >  7.(7), 


over  the  Cr  region,  we  need  a  physical  model  which  gives  us 
a  geometrical  expression  for  the  incremental  volume  dV  and 
for  the  integration  domain.  In  general,  given  a  model  for  the 
region  CR,  a  characteristic  length  Lc  can  be  defined  such  that 
when  7,(7)  >  Lc  oxygen  pervades  the  whole  region  while 
when  7,(7)  <  7,c  this  would  not  happen.  Then  two  diffusion 
regimen  can  be  recognized: 

1 .  for  7,(7)  7,c,  the  oxygen  concentration  profile  is  nearly 

flat  and  we  have  3(7)  =  (co^Cd(7)),’Vcr; 

2.  for  7,(7)  <  Lc,  oxygen  diffusion  is  limited;  the  integral 
must  be  explicitly  calculated  and  depends  on  the 
physical  model. 

Many  mechanistic  models  found  in  the  literature  propose 
a  uniform  Cr  region  and  no  limitations  for  the  oxygen 
diffusion.  Within  these  models,  the  characteristic  length  is 
the  Cr  thickness,  Lqr. 

A  more  realistic  model  of  the  reactive  region,  Cr,  results 
if  we  consider  gas  pores  going  through  the  region,  and  both 
longitudinal  (ID)  gas-phase  and  transversal  (2D)  liquid 
phase  oxygen  diffusion  mechanisms.  For  this  model,  we 
have  two  characteristic  lengths:  the  mean  pore  length  ( Lp ) 
which  characterizes  the  longitudinal  diffusion  and  the  mean 
pore-pore  distance  ( Dp )  which  characterizes  the  transversal 
diffusion. 
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Table  1 

Oxygen  concentration  integrals  for  different  diffusion  mechanisms 


Model 


Limited  diffusion  mechanisms 


Conditions 


Integral 


Nd 


Homogeneous 

L»Lc  E 

L<LCr 

Pore  gas 

None 

U  »Lp 
U  »DP 

Longitudinal 

L  <  Lp 

Transversal 

L  <  Dp 
L>  Rp 

L  <  Dp 
Rp 

Both 

U  <  Dp 
U  >  RP 
L\  '  Lp 

U  <  Dp 

L,<^Rp 

L\  '  Lp 

{ccoRJCDwy 


-Cr/Cd 

co2 


(/) 


(<§/CDW)r 


(4:/Cdw)' 


Cr/Cd  / 

fo^)(cCR/cD(/))y 


I 

Cr/Cd 

o2 


(I) 


(c%,Cd(I)Y 


Cr/Cd 

co2 


(/) 


(c£/CDW)y 


~S/Cd(/) 

/ 


(cS/CDW)r 


0 

1 

0 

1 

2 

1 

5 


Cr/Cd 

o2 


(/) 


(c&/Cd(/))j 


3 


In  the  Appendix  A,  the  value  of  3(1)  is  estimated  for  the 
homogeneous  and  gas  pores  models.  Within  the  latter  model 
the  cases  where  only  longitudinal,  only  transversal  or  both 
diffusion  mechanisms  are  limited  are  taken  into  account. 
The  results  are  summarized  in  Table  1.  It  is  seen  that  in  all 
the  considered  cases,  the  value  of  3(1)  can  be  expressed  as 

/  Cr/c d(t\\nM 

s(/)«(4“/CD(/)ri^j  ,  ns) 

which  is  the  searched  equation. 

The  values  of  the  parameter  Na(I),  give  indication  on  the 
oxygen  diffusion  mechanism  which  is  set  on  at  a  given  cell 
current  density.  Note  that  most  of  the  cases  presented  in 
Table  1  (specially  the  cases  with  the  conditions  <Cand//>), 
correspond  to  border-line  regimens.  At  intermediate  condi¬ 
tions,  transition  regimens  join  smoothly  the  border-line 
ones. 

3.2.  Interface  oxygen  concentration 

To  express  the  oxygen  concentration  at  the  Cr/Cd  inter¬ 
face  as  a  function  of  the  current  density,  we  model  the  gas 
diffusion  through  the  pores  of  the  cathode  diffusion  region 
with  the  Stefan-Maxwell  equations 

V.r,  =  ^2  -fi-  (xtNj  -  xjNi)  i  =  N2,  w,  02,  (16) 

l=i  Uj-i 

where  Dffj  is  an  effective  i  —  j  pair  diffusion  coefficient  in 
the  porous  medium,  A,  the  molar  flux  of  species  i.  and  x,  its 
molar  fraction. 


By  following  the  derivation  of  [8],  which  is  based  on  the 
assumptions  that  both  the  flux  of  nitrogen,  /VN,,  and  the 
gradient  of  water  concentration,  V.rw,  are  zero,  Eq.  (16) 
gives  for  the  gradient  of  nitrogen  concentration  the  expres¬ 
sion 


rieff  ^  v  /  jc-  II  I  v  /)c(T 
^N2-02  XOlUvi-^1  ^  xN2-C'w-02 

(17) 

By  knowing  that  /.)/"  «  /J/i1  [17],  that  the  02  flux  is 

related  to  the  cell  current,  /Vo2)g  =I/(seF),  and  that  the 
effective  diffusion  depends  on  the  gas-phase  volume  frac¬ 
tion,  Dffj  =  If-j(f  [9],  Eq.  (17)  can  be  written  as 

A 

VxN2=xn  ,— /,  (18) 

£s 

where  the  constant  term  A  is 


P 

—  Vxn2  =  *N2No2,g 


,PDn2- o2  (1  —  wO^Av-nJ 

In  these  relations  ,vc  is  the  number  of  electrons  exchanged  in 
the  electrochemical  reactions,  u  is  an  empirical  constant  [9], 
and  Di-j  the  pair  diffusion  coefficients  in  a  non-porous 
medium. 

When  £g  is  "-independent,  the  differential  equation  (18) 
can  be  solved  to  give 


x 


Cr/Cd 

n2 


=  *N2  exP 


(20) 
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where  L  is  the  thickness  of  the  diffusive  region  and  the  zero 
of  the  "-coordinate  is  fixed  at  the  Cf/Cd  interface. 

From  the  relation  xq2  +  xN,  =  1  —  xw,  which  is  constant 
through  the  gas  pores  of  the  cathode,  Eq.  (20)  can  be 
rewritten  for  Xo^Cd  as 


Cr/Cd 

■*o2 


(21) 


From  this  relation  we  see  that  at  I  =  0,  x'o^Cd  =  Xq2,  and 
that  the  limiting  current  I\,  i.e.  the  current  density  value  at 
which  the  oxygen  concentration  x£/CD=0is 


I\  =  In 


4J  al 


(22) 


where  we  made  the  important  step  of  having  the  gas  porosity 
eg  as  a  function  of  the  current  density  7. 

Finally,  and  by  knowing  that  at  the  Cr/Cd  interface  the 
dissolved-oxygen  concentration  in  the  membrane  phase  is 
related  to  the  gas-phase  molar  fraction  by  the  Henry’s  law 
constant 


^  Cr/Cd  _  Cr/Cd 

p  C02  A02  5 

we  can  write  for  c^JC°(I)  the  expression 


(23) 


Cr/Cd 

C02 


(I)=c 


Cr/Cd 

02 


(0) 


1- 


0 

'lN2 

V02 


///i(fg(/i)/fg(/)) 


(24) 


For  small  inlet  oxygen-nitrogen  mole  ratio,  first-order 
expansion  gives  the  simpler  expression 


Cr/Cd 

c02 


(/)  =  c 


Cr/Cd 

o2 


(0) 


(25) 


To  characterize  the  concentration  functions  (24)  and  (25), 
we  analyze  two  elementary  functions  for  £g(7). 


1.  The  exponential  function 


eg7) 


£ 


l5(l-///i) 


II 

OQ  O 

for 

7  =  0, 

=  4 

for 

I  =  h, 

(26) 

=  0 

for 

I  =  oo. 

2.  The  linear  function 


e°g  for  7  =  0, 
4  for  /  =  /]. 


(27) 


In  both  cases,  £0  is  a  physical  parameter  which  depends  on 
the  preparation  of  the  diffusion  part  of  the  MEA,  while  e 
and  I\  are  parameters  which  characterize  the  flooding  pro¬ 
cess.  Eq.  (22)  links  £g  with  I\,  therefore,  we  have  just  one 
independent  parameter. 

To  study  the  behavior  of  Eq.  (24)  with  £g  (7)  given  by  the 
exponential  and  lineal'  function,  we  fix  I\  1 , 
co2^Cd(0)  =  1,  H  =  1-5  from  [9]  and  the  inlet  oxygen  nitro¬ 
gen  mole  ratio  to  0.27  which  corresponds  to  air  feeding,  and 
variate  the  parameter  S. 

In  Fig.  3A,  exponential  and  linear  porosity  functions  are 
inserted  into  Eq.  (24)  and  compared  for  two  different  values 
of  S.  In  addition,  the  exponential  porosity  function  is 
inserted  into  the  approximate  Eq.  (25)  and  plotted  as  well. 

The  case  S  =  1  (elg  =  ([])  corresponds  to  the  case  where 
there  is  no  flooding;  in  this  case,  the  two  curves  coming  from 
expression  (24)  coincide  and  are  nearly  linear.  From  the 
figure,  we  see  that  the  qualitative  behavior  of  the  three 
functions  while  changing  S  is  the  same.  This  is  even  more 
evident  when  we  look  at  the  logarithms  of  the  same  func¬ 
tions  (Fig.  3B)  which  is  the  function  we  need  in  our  relations 
(see  Eq.  (10)).  Notice,  however,  that  the  choice  of  the 
porosity  function  brings  slight  changes  in  the  value  of  S. 
In  conclusion,  we  can  select  arbitrarily  one  of  the  three 
expressions  but  bearing  in  mind  that  the  value  of  S  is  not 
strictly  quantitative.  By  using  a  simplicity  criterion,  we 


Fig.  3.  Oxygen  concentration  with  different  porosity  functions. 
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choose  the  exponential  porosity  function,  and  Eq.  (25) 
becomes 


cR/cD 

co2 


(/)  =  cg/CD(  0) 


1 

h 


(28) 


meaning;  on  the  other  hand,  through  the  experimental  values 
of  the  empirical  parameters,  we  obtain  indications  on  the 
values  of  the  parameters  N&  and  S,  which  characterize  the 
physical  diffusion  mechanisms  inside  the  reactive  and  dif¬ 
fusive  regions,  respectively. 


3.3.  The  semi-empirical  equation 


4.1.  Analysis  of  empirical  parameters 


The  results  of  the  previous  subsections  are  now  inserted 
into  Eq.  (10).  By  first  using  the  results  of  Eq.  (15),  we  get 

1 


77  i  nt  =  K~ 


(1+tfdW)  ln  GO 


-(y  +  W)  ln(c£/Co(/)) 


(29) 


The  first  term  of  this  equation  is  dominant  at  low  current 
densities  and  gives  an  expression  for  the  activation  over¬ 
potential,  while  the  second  term  becomes  dominant  where 
the  oxygen  concentration  is  small  and  expresses  the  con¬ 
centration  overpotential.  It  follows  that  when  this  equation  is 
used  as  a  fitting  equation,  the  coefficient  value  of  the  first 
logarithm  term  in  the  square  brackets  will  be  extracted  from 
the  experimental  behavior  at  low  I ,  while  that  of  the  second 
logarithm  term  from  the  experimental  behavior  at  high  1.  For 
this  reason,  we  reduce  the  continuous  I  dependence  of  /V(|  to 
two  limiting  values,  one  at  low  I,  N®,  and  one  at  high  I,  N\, 
and  write  rjfnt  as 


>L  =k- 


Otc^f 


(1+a£)  in  (/)  —  (y  +  tvJ) 


In  ^ 


(30) 


in  which  we  also  made  use  of  expression  (28). 
Finally,  the  cell  potential,  Eq.  (5),  becomes 


Fee  1 1 


Eq  —  Rceiil  —  b  In  (/)  +  a  In  (  1 


with 

l  +  n9 

E0  =  E(I  =  0)  +  K-  b  = - -*■, 

ctcpp 


1 

h  ) 


(31) 


y  +  *fl 

(32) 


which  is  the  main  result  of  this  section.  In  Eq.  (31),  all  the 
/-dependent  terms  have  a  well  precise  origin,  and  all  para¬ 
meters  have  a  specified  physical  meaning. 


4.  Comparisons  between  theory  and  experiments 

The  interest  of  a  comparison  between  Eq.  (31)  and 
empirical  equations  as  that  of  Kim  et  al.  [1]  or  Squadrito 
et  al.  [3]  is  double.  On  one  hand,  we  gain  a  better  insight 
on  the  empirical  parameters  by  giving  them  a  physical 


Let  now  compare  our  expression  for  the  cell  potential, 
Eq.  (31),  with  the  empirical  equation  given  by  Squadrito 
et  al.  [3]  and  Kim  et  al.  [1]: 

Vcen=E0—RcMI—b  In  (I)+a  In  ^1  —  (this  work) 

Kell  =Eo,s  — /feeu/— In  (/)  +asIk  In  ( 1  — /?/)  (Squadrito) 
bv 

Keii=£o,K-/?ceii/-^— In (/) -m exp(n/)  (Kim) 

(33) 

The  equation  of  Squadrito  et  al.  [3]  fits  quite  well  our 
equation.  The  only  apparent  difference  is  in  the  shape  of 
the  concentration  overpotential  (last  logarithm  term  on  the 
right  side). 

To  find  correspondences,  in  Fig.  4,  we  plot  the  argument 
of  the  second  logarithm  term  of  the  Squadrito  relation 

/(/)  =  (1  -  j?/)as/\  (34) 

for  /j  =  1  / / 1  =  1 ,  «s  =  (0.85 //i)*,  and  different  values  of  k. 
For  such  parameter  values,  we  notice  that  the  functional 
form  of  the  curves  in  this  figure  compares  very  well  with  that 
of  the  curves  of  Fig.  3A.  This  similarity,  leads  us  to  interpret 
the  k  parameter  in  terms  of  the  water  flooding  phenomena: 
k  =  0  corresponds  to  S  =  1  which  means  the  absence  of 
water  flooding,  while  k  >  0  corresponds  to  S  >  1  (i.e. 
4  <  which  in  turn  means  that  water  flooding  reduced 
the  gas  porosity.  The  similarity  also  allows  us  to  define  a 
rescaled  quantity  a's  =  (I\/0.S5)kas  which  now  corresponds 
to  our  a  parameter. 


Fig.  4.  /(/)  at  different  values  of  the  k  parameter. 
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Table  2 

Values  of  the  parameters  given  in  Eq.  (32)  and  needed  to  fit  the  Squadrito  et  al.  [3]  experimental  data 


Membrane 

Temperature  (°C) 

Gas  pressures  (bar) 

E0  (V) 

Rceii  (O/cm2) 

/,  (A/cm2) 

is/ In  10  (V) 

a's  (V) 

a's  In  10/ Z?s 

NF117 

70 

l/t 

0.937 

0.30 

0.824 

0.024 

0.131 

5.4 

1.5/1. 5 

0.936 

0.28 

0.880 

0.026 

0.144 

5.5 

2.5/3 

0.962 

0.28 

0.907 

0.025 

0.152 

6.1 

3/5 

0.970 

0.24 

0.962 

0.026 

0.179 

6.9 

NF115 

70 

1/1 

0.933 

0.24 

0.810 

0.025 

0.099 

4.0 

1.5/1 .5 

0.940 

0.23 

0.833 

0.026 

0.102 

3.9 

2.5/3 

0.956 

0.25 

0.923 

0.023 

0.080 

3.5 

NF112 

70 

2.5/3 

0.946 

0.19 

1.006 

0.022 

0.062 

2.8 

NF117 

80 

3/5 

0.961 

0.28 

0.965 

0.023 

0.112 

4.9 

NF115 

0.966 

0.23 

0.972 

0.022 

0.077 

3.5 

NF112 

0.955 

0.18 

1.008 

0.021 

0.070 

3.3 

Other  links  between  Squadrito  et  al.  [3]  parameters  and 
those  appearing  in  our  empirical  equation  can  be  found  by 
comparing  the  terms  with  the  same  /-dependence  type: 


£o,s  =  E{I  =  0)  +  K: 

k  =  S ;  P  =  \. 

h 


bs 

In  10 


1 

«c/?F 


(35) 


From  the  analysis  done  in  the  previous  sections,  it  follows 
that  k  and  fi  depend  on  the  behavior  of  the  oxygen  con¬ 
centration  in  the  Cq  region  while  hs  and  a's  depend  on  the 
diffusion  mechanisms  inside  the  Cr  region. 

Kim  et  al.  LI]  use  the  equation  shown  on  the  third  sub¬ 
equation  of  Eq.  (33)to  characterize  both  cells  with  air  and 
with  pure  oxygen  feeding.  Since  our  analysis  is  for  cells 
with  air  feeding,  it  is  not  possible  to  justify  the  use,  in  this 
case,  of  their  exponential  term.  The  Kim  exponential 
term  could  still  be  the  correct  term  to  describe  the  concen¬ 
tration  overpotential  for  cells  with  pure  oxygen  feeding. 
However,  to  establish  this  point  our  analysis  requires  further 
extensions. 


4.2.  Comparisons  with  experimental  data 

We  start  this  subsection  by  first  analyzing  the  values  of  the 
fitting  parameters  determined  by  Squadrito  et  al.  [3],  in  order 
to  assess  the  phenomenology  which  is  set  on  in  the  cell  as  the 
current  density  is  increased. 

Table  2  reports  the  parameters  obtained  by  fitting  a 
number  of  experimental  polarization  curves  in  different 
conditions  of  hydrogen/air  pressure,  membrane  thickness 
and  temperature.  The  table  shows  the  same  data  as  in  Table  3 
of  [3],  with  the  values  of  the  parameters  b  and  a  extracted 
from  those  of  bs  and  a s  (see  Eqs.  (32)  and  (35). 

Squadrito  et  al.  found  for  the  k  parameter  a  range  of  values 
between  2  and  4  in  the  90%  of  fitted  cases  (for  the  fitting 
curves  of  Table  2,  k  is  fixed  to  3).  This  means  that  flooding 
phenomena  are  present,  and  that  they  reduce  the  porosity  of 
the  Cd  region  as  the  cell  current  density  is  increased. 


To  start  to  see  if  our  semi-empirical  analysis  is  consis¬ 
tent,  we  can  use  Eq.  (22)  to  find  the  relation  between  S 
and  I\ 


i  i\ 

x^ALh) 


t 


(36) 


The  value  of  A,  see  Eq.  (19),  depends  smoothly  on  tem¬ 
perature  and  at  70  °C  is  worth  0.37  cm/A;  the  rate  between 
nitrogen  and  oxygen  molar  fraction  in  air  is  3.76;  the  values 
of  and  L  depend  on  the  particular  preparation  of  the 
diffusion  region,  and  Bernardi  and  Verbrugge  [9]  report 

=  0.4,  L  =  0.026  cm,  and  they  also  give  p  =  1.5.  With 
these  values  and  for  I\  =  0.9  A/cm  (see  Table  2),  we  get 
S  =  5.3. 

In  Fig.  5,  we  compare  the  concentration  profiles  given  by 
Eq.  (34)  for  k  =  2  and  4  with  that  given,  for  S  =  5.3,  by 
Eq.  (28)  and  by  the  corresponding  equation  with  the  linear 
gas  porosity,  Eq.  (27).  We  observe  that  both  curves  are 
consistent  with  the  range  of  k,  confirming  the  validity  of  the 
analysis  about  the  physical  meaning  of  the  k  parameter,  and 
indicating  that,  at  the  limiting  current,  the  gas  porosity  is 
reduced  by  about  a  factor  of  5  with  respect  to  its  value  at 
zero-current  density. 


Table  3 

Fitting  parameter  for  Eq.  (31) 


Gas 

pressures 

(bar) 

02  (%) 

E0  (V) 

^cell 

(Q/cm2) 

h 

(A/cm2) 

b  (V) 

S 

1/1 

27 

0.933 

0.21 

0.816 

0.026 

3.0 

5.7 

2.5/3 

0.956 

0.225 

0.928 

0.023 

2.5 

5.2 

1/1 

0.940 

0.26 

0.763 

0.027 

7.5 

5.9 

1/3 

0.982 

0.26 

1.070 

0.028 

5.5 

4.7 

1/5 

0.996 

0.29 

1.149 

0.027 

5.5 

4.5 

1/1 

20 

0.917 

0.38 

0.789 

0.020 

4.5 

5.5 

5 

0.879 

0.88 

0.309 

0.020 

1.5 

3.6 
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The  last  column  of  Table  2  shows  the  rates  between  the  a 
and  the  b  parameters.  From  Eq.  (32),  we  have 


a  _  1  + 

b  1  +  Nl  ‘ 


(37) 


Because  the  values  found  in  literature  for  y  are  in  the  range 
0.5-1,  we  get  agreement  with  the  experimental  data  for 
IVj  =  0  and  Nld  in  the  range  2-5.  As  it  is  clear  from  Eq.  (13), 
at  very  low  I  the  diffusion  oxygen  length  becomes  very 
large,  and  the  reactive  region  is  fully  permeated  by  the 
oxygen  molecules,  VCff  =  FcR-  The  choice  of  the  fitting 
value  Nd  =  0  is  then  perfectly  justified  (see  Table  1,  “none” 
case). 

Instead,  the  range  of  Nld  values  tells  us  that  at  high  I  the 
diffusion  oxygen  length  becomes  smaller  than  the  charac¬ 
teristic  lengths  of  the  Cr  porous  region,  and  the  oxygen 
molecules  are  unable  to  permeate  all  VcE:  concentration 
overpotetial  is  set  on.  Once  more,  the  fitting  values  of  are 
justified  by  our  model.  We  then  believe  that  the  diffusion 
oxygen  phenomenology,  on  which  our  model  as  been  con¬ 
structed,  is  substantially  correct. 

The  last  information  we  can  extract  from  Table  2  is  on  the 
value  of  the  cathodic  transfer  coefficient  ac.  Indeed,  once  the 
diffusion  model  is  defined,  the  values  of  Nd  are  known,  and 
from  the  third  sub-equation  of  Eq.  (35),  we  can  determine 


the  value  of  the  cathodic  transfer  coefficient  ac.  A  value  of 
Nd  —  0  would  mean  ac  in  the  range  1.1-1. 4  which  seems  to 
indicate  a  two  electron  rate  determining  step  for  the  cathodic 
reaction. 

Finally,  we  use  our  Eq.  (31)  as  an  ordinary  semi-empirical 
equation  to  fit  experimental  polarization  curves.  The  value 
of  Nfj  has  been  fixed  to  0  and  S  has  been  linked  to  I\  through 
Eq.  (36). 

In  Fig.  6,  we  show  the  comparison  between  experimental 
data  and  fitted  polarization  curves.  The  experimental  data 
of  Fig.  6A  are  generated  by  inserting  the  fitting  parameters 
reported  in  Table  3  of  Squadrito  et  al.  [3]  into  the  corre¬ 
sponding  empirical  equation,  which  is  given  in  the  second 
sub-equation  of  Eq.  (33).  The  experimental  data  of  Fig.  6B 
and  C  are  obtained  through  the  same  procedure  but  using 
the  fitting  data  of  Tables  2  and  3  of  Kim  et  al.  [1]  and 
the  fitting  equation  in  the  third  sub-equation  of  Eq.  (33). 
The  experiments  of  Fig.  6A  and  B  are  performed  at  70  °C 
and  with  a  Nation  H  117  membrane,  while  the  experiments 
of  Fig.  6C  are  performed  at  60  °C  and  with  an  Asahi  Che¬ 
mical  Aciplex-S  membrane.  A  good  agreement  between 
experimental  data  and  our  fitting  results  is  found  in  all  three 
cases. 

Table  3  gives  the  values  of  the  fitting  parameters.  The 
value  of  S  is  shown  even  if  is  not  an  independent  parameter 
because  of  its  physical  meaning.  The  behaviour  of  S  fits  very 
well  with  what  we  expect  about  the  water  flooding  phenom¬ 
ena.  Indeed,  by  increasing  the  gas  pressure,  from  simple 
arguments  of  capillarity,  we  expect  that  water  would  find  a 
larger  resistance  to  flooding  and  we  observe,  as  expected, 
a  decline  of  S  [22],  Furthermore,  the  net  decrease  of  S 
observed  when  the  oxygen  percentage  drops  from  20  to 
5,  is  easily  related  to  the  corresponding  fall  of  the  limiting 
current:  by  decreasing  the  current  density  the  water  gener¬ 
ated  in  the  cathode  reactive  region  decreases  as  well;  this 
generates  a  pressure  drop,  and  the  flooding  is  also  expected 
to  decrease. 

The  observed  fall  of  N\  when  the  oxygen  percentage 
drops  from  20  to  5,  is  a  very  encouraging  indication  of  the 
validity  of  its  physical  interpretation.  Indeed,  by  following 
the  same  arguments  expressed  after  Eq.  (29),  it  appears 
evident  that  when  I\  approaches  0,  also  Nld  must  approach 
Nd,  and  this  is  clearly  confirmed  by  the  results. 


0.0  0.2  0.4  0.6  0.8  1.0  0.0  0.2  0.4  0.6  0.8  1.0  0.0  0.2  0.4  0.6  0.8 


Fig.  6.  Comparison  between  experimental  data  and  theoretical  results.  Experimental  data  are  from  [3]  (A),  and  from  [1]  (B  and  C). 
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However,  the  values  of  /Vd  appear  to  be  slightly  too  high 
for  the  curves  represented  in  Fig.  6B.  There  could  be  two 
reasons  to  explain  this  behavior: 

1.  We  are  not  fitting  directly  experimental  data  but  values 
which  in  turn  are  obtained  by  a  fit  of  the  experiments. 
It  is  quite  likely  that  a  high  value  of  Nd  is  needed 
to  represent  the  exponential  fall  which  appears  in 
the  fitting  expression  (third  sub-equation  in  Eq.  (33)) 
used  to  represent  the  experimental  data  of  cases  B  and  C 
but  not  necessarily  on  the  fall  of  the  true  experimental 
data. 

2.  We  are  considering  the  lack  of  oxygen  due  to  the  limited 
oxygen-in-nitrogen  diffusion  rate,  as  the  only  cause  for 
the  non-linear  drop  of  the  polarization  curve  at  high 
current  density.  Fig.  10  and  Table  3  of  Kim  et  al.  [1] 
(corresponding  to  case  C)  confirms  this  hypothesis; 
indeed  the  data  obtained  by  using  pure  oxygen,  and 
shown  in  the  same  figure,  do  not  present  any  trace  of 
non-linear  drop.  On  the  other  hand.  Fig.  5  and  Table  2  of 
the  same  paper  (from  which  we  took  the  air  case  B) 
show  a  consistent  drop  even  in  the  case  of  pure  oxygen. 
It  is  evident  that  for  this  last  preparation  other  causes  of 
non-linearity  (as,  e.g.  non-linear  Ohmic  phenomena) 
have  to  be  taken  into  account.  To  verify  this  hypothesis, 
we  added  the  non-linear  fall  of  the  pure  oxygen  to  the 
fitting  curve  for  air:  we  observe  that,  by  this  procedure, 
Nd  drops  to  a  value  around  4  (well  inside  the  expected 
range)  in  the  case  of  1/3  and  1/5  atm  pressures,  while  it 
stays  very  high  in  the  case  of  1/1  atm  pressures. 

4.3.  Comparisons  with  published  mechanistic  models 

In  this  new  light,  it  is  also  relevant  to  review  other  models 
found  in  the  literature.  Because  all  of  the  model  examined 
suppose  an  uniform  reactive  region,  the  only  characteristic 
length  is  its  thickness  and  we  may  have  for  Nd  values  of  0  or 
1  (see  Table  1  homogeneous  case).  Some  of  the  authors 
(Amphlett  et  al.  [4,5],  Mann  et  al.  [6],  Kulikovsky  et  al. 
[12,13],  Garau  [14])  assume  oxygen  diffusion  in  gas-phase; 
in  this  case,  even  at  high  current  density,  oxygen  reaches  the 
whole  reactive  region  and  we  have  that  both  and  Nd  are 
worth  0.  Other  authors  (Bernardi  and  Verbrugge  [9],  Rho 
et  al.  [18])  assume  oxygen  diffusion  in  liquid  phase;  this 
imply  that  already  at  low  current  density,  oxygen  did  not 
reaches  the  whole  reactive  region  and  we  have  that  both  Nd 
and  Nd  are  worth  1.  Some  other  authors  (Springer  et  al.  [10], 
Marr  and  Li  [19],  Um  et  al.  [20])  use  intermediate  values  for 
oxygen  diffusion  coefficients;  in  these  cases  we  may  have 
that  Nd  =  0  and  Nd  =  1 .  In  all  of  these  papers,  the  value  of  y 
is  inside  the  range  0.5-1. 

These  values  of  /V(|  and  y  imply  a/b  ratios  in  the  range 
0.5-2,  which  is  considerably  lower  than  the  range  shown  in 
Table  2.  This  is  the  reason  why  the  values  and  the  shape  of 
the  performance  curves  given  by  these  models,  near  the 


limiting  current,  are  quantitatively  wrong  and  qualitatively 
different  from  the  experimental  ones. 

Only  Springer  et  al.  [10]  tried  to  solve  this  discrepancy, 
and  to  reach  this  goal  they  used  the  ionic  conductivity  and 
the  oxygen  diffusion  coefficient  in  the  reactive  region  as 
fitting  parameters.  They  observed  that  by  using  very  low 
conductivity  values,  and  ad  hoc  values  for  the  diffusion 
coefficient,  the  modelled  curves  get  closer  to  the  experi¬ 
mental  ones.  This  is  due  to  the  fact  that  when  the  current 
density  increases,  oxygen  tends  to  concentrate  near  the 
interface  with  the  gas  diffuser,  then  protons  have  to  travel 
a  longer  way  through  the  reactive  region,  and  the  Ohmic 
losses  increase  more  than  linearly.  This  non-linear  contri¬ 
bution  (neglected  by  us  due  to  the  small  thickness  of  the 
reactive  region),  becomes  important  when  the  proton  con¬ 
ductivity  gets  very  low.  Springer  et  al.  [10],  justified  the  very 
low  values  of  proton  conductivity  (^0.001  S/cm)  by  experi¬ 
mental  measurements  made  at  low  humidification  degree. 
However,  it  is  very  hard  to  believe  that  such  conditions 
could  exist  in  the  cathode  reactive  region,  where  water  is 
generated  by  the  reaction. 

We  believe  that  the  correct  way  to  get  agreement  with  the 
experimental  data,  is  to  change  the  diffusion  model  for  the 
oxygen  in  the  reactive  cathode  region  in  such  a  way  that 
=  0  and  Nld  €  1-5.  This  can  be  achieved  by  supposing 
the  existence  of  a  small  amount  of  gas  pores  in  the  reactive 
cathode  region,  and  by  taking  into  account  both  for  the 
longitudinal  (gas-phase)  oxygen  diffusion  along  the  pores 
and  for  the  transversal  (liquid  phase)  oxygen  diffusion  inside 
the  membrane  phase  [21]. 

5.  Summary  and  conclusions 

The  presented  semi-empirical  derivation  of  a  perfor¬ 
mance  equation  for  PEMFCs  is  based  on  the  observation 
that  the  strongest  non-linear  contributions  to  the  cell  poten¬ 
tial  drop  at  high  current  density  are  arising  from  interface 
phenomena  happening  in  the  cathode  reactive  region.  The 
key  quantity  to  describe  such  contributions  appears  to  be 
the  integral  of  the  oxygen  concentration  over  the  reactive 
region.  To  evaluate  such  integral,  oxygen  propagation  from 
the  inlet  gas  channel  trough  the  diffusion  region  of  the 
cathode,  and  then  inside  the  reactive  region,  has  been 
described. 

The  oxygen  propagation  inside  the  diffusion  region  is 
determined  by  integrating  the  Stefan-Maxwell  equations, 
with  the  assumption  that  both  the  nitrogen  flux  and  the  water 
concentration  gradient  are  zero.  We  also  assume  that  the  gas 
porosity  of  the  region  is  a  decreasing  function  of  the  cell 
current  density,  in  a  way  to  simulate  the  well  known  flooding 
phenomenon  of  the  cathode.  The  shape  of  this  function  has 
been  found  to  be  of  minor  importance  on  the  characteriza¬ 
tion  of  flooding,  while  the  important  parameter  is  recognized 
to  be  the  ratio  between  the  gas  porosity  at  zero  and  at 
limiting  current  densities. 
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Permeation  of  oxygen  molecules  inside  the  reactive 
region  can  be  qualitatively  described  by  comparing  the 
distance  travelled  by  the  oxygen  molecules  with  the  char¬ 
acteristic  length  scales  of  the  region.  If  the  travelled  distance 
is  larger  than  the  region  length  scales  (such  as  the  mean  pore 
length  or  the  mean  pore  distance)  oxygen  permeates  the 
whole  region.  The  more  the  travelled  distance  reduces,  the 
more  the  oxygen  permeated  volume  depends  on  the  details 
of  the  geometrical  structure  of  the  reactive  region,  and  the 
value  of  the  integral  of  the  oxygen  concentration  changes, 
too.  The  quantitative  evaluation  of  such  integral  for  two 
different  models  and  for  different  ranges  of  the  travelled 
distance,  brings  to  an  expression  containing  a  parameter  Na 
which  characterizes  the  structure  of  the  reactive  region. 

The  final  empirical  equation  for  the  performance  curve 
looks  similar  to  other  published  equations  (see,  e.g.  Squa- 
drito  et  al.  [3]),  but  in  our  case  all  equation  parameters  have  a 
clear  physical  origin  and  a  clear  physical  meaning. 

The  empirical  equation  is  then  used  to  fit  a  number  of 
experimental  performance  curves.  All  analyzed  experimen¬ 
tal  data  are  well  reproduced,  especially  at  high  current 
density. 

The  model  used  to  derive  the  empirical  equation,  and  the 
set  of  numerical  values  found  for  the  fitting  parameters  give 
strong  evidences  of  the  existence  of  flooding  phenomena 
inside  the  cathode  diffusion  region,  and  of  the  presence  of  a 
gas  pores  structure  inside  the  cathode  reactive  region. 

Then,  the  modelling  of  such  phenomena  and  structures 
seems  to  be  necessary  in  order  to  get  a  quantitative  descrip¬ 
tion  of  the  performance  drop  at  high  current  density  of  an 
Ha/air  feed  cell.  This  kind  of  modelling  is  being  used  by  us 
within  a  mechanistic  model  based  on  a  Bernardi  and  Ver- 
brugge  type  of  approach  |7,9].  Work  in  this  direction  is  in 
progress  [21]. 
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Appendix  A.  Evaluation  of  oxygen 
concentration  integrals 

From  the  equation  for  the  oxygen  concentration  profile, 
Eq.  (14),  the  integral  in  Eq.  (11)  is  written  as 

m  =  (co:/Cd(/))?  Jc  (i  -^) '  dv.  (A.i) 


•  Homogeneous  model 

In  this  case,  we  have  dV  =  gq  dZ  and  the  integral 
becomes 

5(7)  =  L(/)c7C(c£/Cd(/))1’  f  (1  -x)y  d v.  (A. 2) 

Jo 

From  Eq.  (13)  and  because  in  this  case  a  =  a c  is  a 
constant,  we  get 

Cr/CD/i\ 

S(I)cx^j^(cC0fD(I)Y,  (A. 3) 

which  is  the  expression  reported  in  Table  1,  second  row. 

•  Pore  gas  model 

o  Longitudinal  limitation 

In  this  case,  we  have  dV  =  VcR  dl\/ Lp  and,  by  fol¬ 
lowing  the  same  steps  as  for  the  homogeneous  model, 
we  get 

CR/CD/r\ 

S(/)  a (/))»,  (A. 4) 

which  is  the  expression  reported  in  Table  1,  fourth  row. 
o  Transversal  limitation 

In  this  case,  we  have  dV  =  2nLp(lt  +  Rp)  d/t,  where 
Rp  is  the  mean  pore  radius,  and  the  integral  becomes 

5(7)  =  27tLfLp(c^/CD(7))y  J  (1  -  x)y  (x  +  ^  dx. 

(A. 5) 

The  diffusion  section  a  is  given  by  the  internal  surface 
of  the  pores,  er  =  2L(  R< / Rp,  which,  once  more,  is  a 
constant. 

In  the  case,  where  Lt^>  Rp,  the  ratio  Rp/Lt  inside  the 
integral  can  be  neglected  and,  we  get 

5(7)  a  ( cC0fD(I))\  (A. 6) 

which  is  the  expresssion  reported  in  Table  1,  fifth  row. 
On  the  other  extreme,  when  Lt  <C  Rp,  we  can  write 

5(7)  «  27rLtLp(cg/CD(7))'’7?p  [  (1  -  x)y  dx 

Jo 

Cr/Cd/.x 

oc  02  7  U(cg/Cp(7))y,  (A. 7) 

which  is  the  expression  reported  in  Table  1,  sixth  row. 
o  Both  limitations 

In  this  case,  we  have  dV  =  2n(lt  +  Rp)  d/t  d/j  and  the 
integral  becomes 

5(7)  =271 J  dl\Ly(co2(l\:I))y J  (l-x)?(^  +  x)dx. 

(A.  8) 

In  the  case,  where  7.,  Rp,  the  integral  on  x  becomes 
constant  and,  we  have 
rL i 

5(7)  cx  /  d/1Lf(c0,(/i,7))5'. 

Jo 


(A.9) 
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The  value  of  Lt  can  be  estimated  by  using  Eq.  (13)  and 
by  observing  that  a  =  2L\Lcv(snc / Ry,Lv  is  proportional 
to  L\ 

Lt  oc  c02(/i,/)y ,  (A.  10) 

and,  we  have 

S(J)«  (cg/C°(/))y,  (A.  11) 

which  is  the  expression  reported  in  Table  1,  seventh 
row. 

In  the  case  where  Lt  <C  Rp,  as  for  Eq.  (A. 7),  we  get 

rh 

S(7)  cx  /  d/1Lt(Co,(/i,/))7,  (A.  12) 

Jo 

and  finally,  we  have 

m  cx  (cC0fD(I))\  (A.13) 

which  is  the  expression  reported  in  Table  1 ,  eighth  row. 
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